El código y los archivos orginales de esta PEC, así como los informes en formato HTML pueden consultarse en el siguiente repositorio de GitHub. El repositorio es privado, así que se tendrá que requerir el acceso para poder ver el código.
En un estudio sobre una enfermedad del hígado se recogieron datos de 583 pacientes del departamento de digestivo de un hospital. De ellos, 416 experimentan una afección del hígado y 167 tienen otros problemas pero no se ven afectados por dicha afección. El archivo dataset.csv contiene un conjunto de variables que se consideran “marcadores” de dicha enfermedad, además del género de los pacientes. Las variables registradas son:
data <- read.csv("dataset.csv") %>% dplyr::mutate(Selector = factor(Selector))
head(data)
Realizar un resumen numérico y gráfico de los datos. Observar si hay valores faltantes y tenerlo en cuenta en el análisis.
La manera más fácil de realizar un resumen numérico es usar la función summary().
summary(data)
## Age Gender TB DB
## Min. : 4.00 Female:142 Min. : 0.400 Min. : 0.100
## 1st Qu.:33.00 Male :441 1st Qu.: 0.800 1st Qu.: 0.200
## Median :45.00 Median : 1.000 Median : 0.300
## Mean :44.75 Mean : 3.299 Mean : 1.486
## 3rd Qu.:58.00 3rd Qu.: 2.600 3rd Qu.: 1.300
## Max. :90.00 Max. :75.000 Max. :19.700
##
## Alkphos Sgpt Sgot TP
## Min. : 63.0 Min. : 10.00 Min. : 10.0 Min. :2.700
## 1st Qu.: 175.5 1st Qu.: 23.00 1st Qu.: 25.0 1st Qu.:5.800
## Median : 208.0 Median : 35.00 Median : 42.0 Median :6.600
## Mean : 290.6 Mean : 80.71 Mean : 109.9 Mean :6.483
## 3rd Qu.: 298.0 3rd Qu.: 60.50 3rd Qu.: 87.0 3rd Qu.:7.200
## Max. :2110.0 Max. :2000.00 Max. :4929.0 Max. :9.600
##
## ALB A.G Selector
## Min. :0.900 Min. :0.3000 1:416
## 1st Qu.:2.600 1st Qu.:0.7000 2:167
## Median :3.100 Median :0.9300
## Mean :3.142 Mean :0.9471
## 3rd Qu.:3.800 3rd Qu.:1.1000
## Max. :5.500 Max. :2.8000
## NA's :4
El resumen numérico nos puede dar varias informaciones, donde tenemos que la edad de los pacientes oscila entre 4 y 90 años, siendo la media y mediana muy similares (44.75 y 45), entre otras variables.
Buscando si hay valors faltantes (indicados como NA), se puede observar fácilmente que en el ratio de albúmina y globulina (variable A.G), hay un total de 4 valores que no están presentes.
Para proceder con el análisis, dado que el número de valores NA és muy peqeño y solo afecta a una variable, podríamos usar el método de complete-case analysis usando solo los casos que tengan todas las variables completas o el método de available-case analysis. Debido a su mayor senzillez, vamos a proceder con el filtraje de los casos incompletos, quedandonos con esos que tengan información para todas las variables.
data_filt <- data[which(complete.cases(data)),]
summary(data_filt)
Para hacer el resumen gráfico de las variables, podemos usar la función plot(), que nos dibujará una serie gráficos de dispersión relacionando las variables una contra una.
plot(data, main = "Scatterplot comparing all the variables one-vs-one")
Figura 1.1. Gráfico de dispersión que compara todas las variables (por pares) de nuestro set de datos. Los datos usados para este gráfico incluyen los casos incompletos y todas las variables, sean numéricas o no.
A simple vista, las variables indicativas de la bilirubina directa y la bilirubina total (DB y TB, respectivamente) parecen tener una correlación bastante buena, así como las variables TP y ALB o ALB y A/G.
Estudiar la relación entre las variables numéricas, especialmente entre la bilirubina total y la directa. ¿Están correlacionadas?
Para hacerlo sencillo, vamos a crear otro set de datos a partir de los datos filtrados (casos completos) con solo las variables numéricas (todas menos Gender y Selector)
data_num <- data_filt %>% dplyr::select(everything(), -Gender, -Selector)
En el resumen gráfico del apartado anterior, a simple vista parece que la bilirubina total (TB) y la directa (DB) tienen una buena correlación.
Seguidamente vamos a comprobar la correlación entre variables usando la función cor() y a dibujar un gráfico de correlación con la función corrplot().
corrplot::corrplot(corr = cor(data_num), type = "lower",
number.digits = 2, method = "color",
addgrid.col = "gray50", addCoef.col = "black",
title = "Correlation plot showing the correlation between each pair of variables")
Figura 1.2. Correlaciones entre las distintas variables. A más intensidad de color, más correlación, ya sea positiva (azul) o negativa (rojo). Para hacer este gráfico solo se han usado las variables numéricas y los casos completos.
En el gráfico anterior, donde a más intensidad de color más correlación, vemos que, en general, la correlación entre las distintas variables es pequeña (coeficiente de correlación de Pearson de menos de 0.3 en valor absoluto). Sin embargo, hay cuatro pares de variables que sí que tienen una alta correlación, siendo la bilirubina total i la bilirubina directa las variables más correlacionadas (Pearson = 0.87):
Mirando estos números (figura 1.2), además del gráfico de dispersión (figura 1.1), podemos decir que la bilirubina total y la directa están altamente correlacionadas.
Hacer un gráfico multivariante de las correlaciones dos a dos.
Mirar apartado anterior
Obtener una lista de las correlaciones dos a dos, ordenadas de mayor a menor.
En la Tabla 1.1 se muestran las correlaciones dos a dos, ordenadas de mayor a menor valor de Pearson en valor absoluto. Por razones de espacio, la Tabla 1.1 se encuentra cortada a 10 entradas y la lista completa puede verse en la Tabla S1.1 del material suplementario.
cor_list_1d <- cor(data_num) %>%
reshape2::melt() %>%
dplyr::mutate(abs_pearson = sqrt(value^2)) %>%
dplyr::arrange(desc(abs_pearson)) %>%
magrittr::set_colnames(c("Variable 1", "Variable 2", "Pearson's correlation", "Absolute Pearson's"))
cor_list_1d %>% head(10) %>% knitr::kable(align = "c", caption = "Tabla 1.1. Lista de correlaciones dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos. Solo se muestran las 10 primeras entradas. La lista con todos los valores se puede encontrar en el material suplementario (Tabla S1.1).")
| Variable 1 | Variable 2 | Pearson’s correlation | Absolute Pearson’s |
|---|---|---|---|
| Age | Age | 1.000000 | 1.000000 |
| TB | TB | 1.000000 | 1.000000 |
| DB | DB | 1.000000 | 1.000000 |
| Alkphos | Alkphos | 1.000000 | 1.000000 |
| Sgpt | Sgpt | 1.000000 | 1.000000 |
| Sgot | Sgot | 1.000000 | 1.000000 |
| TP | TP | 1.000000 | 1.000000 |
| ALB | ALB | 1.000000 | 1.000000 |
| A.G | A.G | 1.000000 | 1.000000 |
| DB | TB | 0.874481 | 0.874481 |
Repetir el punto anterior considerando las correlaciones parciales entre pares de variables cuando se fijan todas las demás. ¿Se ven afectadas vuestras conclusiones?
Para calcular las correlaciones parciales, usaremos la función pcor() del paquete ppcor().
Primeramente haremos un gráfico de correlación como el de la Figura 1.2, con tal de observar las diferencias. Este gráfico podrá ser observado en la Figura S1.1 en el material suplementario.
par_cor_1e <- ppcor::pcor(data_num)
Viendo la Figura S1.1, se ve claramente que las correlaciones parciales cambian comparadas con las correlaciones obtenidas en la Fig. 1.2. Por ejemplo vemos como la correlación entre la bilirubina total (TB) y la directa (DB) disminuye ligeramente, mientras que las correlaciones parciales ALB vs TP (Pearson’s = 0.82), ALB vs A.G (Pearson’s = 0.90), y TP vs A.G (Pearson’s = -0.69), augmentan (en valor absoluto) mucho comparadas con las correlaciones anteriores.
Igual que en el apartado d, la lista completa de las correlaciones parciales se puede observar en la Tabla S1.2 en el material suplementario mientras que la Tabla 1.2 solo muestra las diez primeras entradas.
par_cor_1e_list <- par_cor_1e$estimate %>%
reshape2::melt() %>%
dplyr::mutate(abs_pearson = sqrt(value^2)) %>%
dplyr::arrange(desc(abs_pearson)) %>%
magrittr::set_colnames(c("Variable 1", "Variable 2", "Pearson's correlation", "Absolute Pearson's"))
par_cor_1e_list %>% head(10) %>% knitr::kable(align = "c", caption = "Tabla 1.2. Lista de correlaciones parciales dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos. Solo se muestran las 10 primeras entradas. La lista con todos los valores se puede encontrar en el material suplementario (Tabla S1.2).")
| Variable 1 | Variable 2 | Pearson’s correlation | Absolute Pearson’s |
|---|---|---|---|
| Age | Age | 1.0000000 | 1.0000000 |
| TB | TB | 1.0000000 | 1.0000000 |
| DB | DB | 1.0000000 | 1.0000000 |
| Alkphos | Alkphos | 1.0000000 | 1.0000000 |
| Sgpt | Sgpt | 1.0000000 | 1.0000000 |
| Sgot | Sgot | 1.0000000 | 1.0000000 |
| TP | TP | 1.0000000 | 1.0000000 |
| ALB | ALB | 1.0000000 | 1.0000000 |
| A.G | A.G | 1.0000000 | 1.0000000 |
| ALB | TP | 0.8971925 | 0.8971925 |
Calcular el vector de medias, la matriz de covarianzas y la matriz de correlaciones de las variables respuesta numéricas conjuntamente pero por separado para los dos grupos (niveles del factor Selector sin tener en cuenta el género).
Primero separamos los datos según el valor en la variable Selector (1 o 2).
data_num_selector1 <- data_filt %>% dplyr::filter(Selector == 1) %>% dplyr::select(everything(), -Gender, -Selector)
data_num_selector2 <- data_filt %>% dplyr::filter(Selector == 2) %>% dplyr::select(everything(), -Gender, -Selector)
mean_vect_selector1 <- colMeans(data_num_selector1)
mean_vect_selector2 <- colMeans(data_num_selector2)
cov_mat_selector1 <- cov(data_num_selector1)
cov_mat_selector2 <- cov(data_num_selector2)
corr_mat_selector1 <- cor(data_num_selector1)
corr_mat_selector2 <- cor(data_num_selector2)
Después de haber separado los datos, calculamos los vectores de medias con colMeans(), la matriz de covarianzas con cov() o la matriz de correlaciones con cor():
Vector de medias:
46.1449275, 4.1804348, 1.9316425, 319.5362319, 99.97343, 138.173913, 6.4586957, 3.0584541, 0.9141787
41.3636364, 1.1448485, 0.3963636, 220.6848485, 33.8363636, 40.7636364, 6.5393939, 3.3393939, 1.0295758
Matriz de covarianzas:
246.1871776, -2.6722181, -1.7319577, 376.866407, -436.8823385, -243.4586799, -3.1780187, -3.4973541, -1.3308976, -2.6722181, 51.2423524, 19.9290225, 318.8957469, 279.9389462, 509.423992, 0.0609333, -1.2364565, -0.4557098, -1.7319577, 19.9290225, 10.3203765, 164.7544198, 137.6199711, 247.9961785, 0.0670728, -0.5629437, -0.1916362, 376.866407, 318.8957469, 164.7544198, 7.227736310^{4}, 5390.7891006, 1.315774210^{4}, -5.9794926, -29.9190722, -18.3494859, -436.8823385, 279.9389462, 137.6199711, 5390.7891006, 4.546146410^{4}, 5.675900710^{4}, -10.6640541, -0.8768935, 1.9592638, -243.4586799, 509.423992, 247.9961785, 1.315774210^{4}, 5.675900710^{4}, 1.143361110^{5}, -7.4678598, -18.6288346, -5.948743, -3.1780187, 0.0609333, 0.0670728, -5.9794926, -10.6640541, -7.4678598, 1.2040283, 0.6576987, 0.0720592, -3.4973541, -1.2364565, -0.5629437, -29.9190722, -0.8768935, -18.6288346, 0.6576987, 0.6200131, 0.170283, -1.3308976, -0.4557098, -0.1916362, -18.3494859, 1.9592638, -5.948743, 0.0720592, 0.170283, 0.1063755
291.0133038, 0.2049335, 0.1732816, -190.6956763, -46.7998891, -61.6452328, -3.2686807, -2.2266075, -0.2056375, 0.2049335, 1.0193178, 0.5197982, 80.9239763, 8.5153104, 14.2240798, -0.1654361, -0.1452531, -0.0472918, 0.1732816, 0.5197982, 0.2724257, 41.6836031, 4.3908647, 7.3869401, -0.079429, -0.0732705, -0.0248308, -190.6956763, 80.9239763, 41.6836031, 2.00301210^{4}, 1341.8261641, 1384.5409091, -4.3997044, -16.1198263, -9.8257443, -46.7998891, 8.5153104, 4.3908647, 1341.8261641, 632.332816, 610.3452328, 0.9814856, 0.8766075, 0.0663326, -61.6452328, 14.2240798, 7.3869401, 1384.5409091, 610.3452328, 1336.8645233, -4.0711197, -2.3125831, 0.2007528, -3.2686807, -0.1654361, -0.079429, -4.3997044, 0.9814856, -4.0711197, 1.1094752, 0.7056338, 0.0987973, -2.2266075, -0.1452531, -0.0732705, -16.1198263, 0.8766075, -2.3125831, 0.7056338, 0.6061826, 0.1649558, -0.2056375, -0.0472918, -0.0248308, -9.8257443, 0.0663326, 0.2007528, 0.0987973, 0.1649558, 0.0825138
Matriz de correlaciones:
1, -0.0237917, -0.0343603, 0.0893416, -0.13059, -0.0458882, -0.1845888, -0.2830782, -0.2600706, -0.0237917, 1, 0.86661, 0.165704, 0.1834117, 0.2104617, 0.0077575, -0.2193632, -0.195188, -0.0343603, 0.86661, 1, 0.1907604, 0.2009148, 0.2282998, 0.0190275, -0.2225444, -0.182898, 0.0893416, 0.165704, 0.1907604, 1, 0.0940437, 0.14474, -0.0202696, -0.141334, -0.2092676, -0.13059, 0.1834117, 0.2009148, 0.0940437, 1, 0.7872658, -0.0455808, -0.0052231, 0.0281741, -0.0458882, 0.2104617, 0.2282998, 0.14474, 0.7872658, 1, -0.0201273, -0.069967, -0.0539402, -0.1845888, 0.0077575, 0.0190275, -0.0202696, -0.0455808, -0.0201273, 1, 0.7612165, 0.2013494, -0.2830782, -0.2193632, -0.2225444, -0.141334, -0.0052231, -0.069967, 0.7612165, 1, 0.6630558, -0.2600706, -0.195188, -0.182898, -0.2092676, 0.0281741, -0.0539402, 0.2013494, 0.6630558, 1
1, 0.0118988, 0.0194613, -0.0789846, -0.1090977, -0.0988324, -0.1819103, -0.167643, -0.0419645, 0.0118988, 1, 0.9864065, 0.5663444, 0.3354075, 0.3853237, -0.1555667, -0.184786, -0.1630677, 0.0194613, 0.9864065, 1, 0.5642862, 0.3345439, 0.3870765, -0.1444763, -0.1803032, -0.1656165, -0.0789846, 0.5663444, 0.5642862, 1, 0.377035, 0.2675595, -0.0295136, -0.1462908, -0.2416909, -0.1090977, 0.3354075, 0.3345439, 0.377035, 1, 0.6638332, 0.0370555, 0.0447745, 0.0091831, -0.0988324, 0.3853237, 0.3870765, 0.2675595, 0.6638332, 1, -0.1057089, -0.0812366, 0.0191141, -0.1819103, -0.1555667, -0.1444763, -0.0295136, 0.0370555, -0.1057089, 1, 0.8604365, 0.3265298, -0.167643, -0.184786, -0.1803032, -0.1462908, 0.0447745, -0.0812366, 0.8604365, 1, 0.7375689, -0.0419645, -0.1630677, -0.1656165, -0.2416909, 0.0091831, 0.0191141, 0.3265298, 0.7375689, 1
Calcular la varianza generalizada, la varianza total y el coeficiente de dependencia global \(\eta^2 = 1−|R|\) para cada grupo por separado.
Varianza generalizada: se defene como el determinante de la matriz de varianzas. Así pues, con la función det() y la función var() vamos a calcular la varianza generalizada (ex. det(var(data))).
var_g_total <- det(var(data_num))
var_g_selector1 <- det(var(data_num_selector1))
var_g_selector2 <- det(var(data_num_selector2))
La varianza generalizada del total de los datos es 5.98808510^{15}, mientras que para los datos etiquetados como Selector == 1 es 3.202195510^{16} y para los datos etiquetados con Selector == 2 es 1.053902710^{7}.
Varianza total: se define como la suma de las varianzas de cada variable o como la traza de la matriz de varianzas y covarianzas. Como las varianzas son la diagonal de la matriz de varianzas-covarianzas, con la suma de la diagonal vamos a poder calcular la varianza total.
vt_global <- sum(diag(var(data_num)))
vt_selector1 <- sum(diag(var(data_num_selector1)))
vt_selector2 <- sum(diag(var(data_num_selector2)))
La varianza total del total de los datos es 1.772031510^{5}, mientras que para los datos etiquetados como Selector == 1 es 2.323846210^{5} y para los datos etiquetados con Selector == 2 es 2.22934210^{4}.
Coeficiente de dependencia global o \(\eta^2 = 1−|R|\): para calcular el coeficiente de dependencia global tenemos que calcular el determinante \(||\) de la matriz de correlaciones \(R\) (\(|R|\)) y restarlo de 1.
coef_dep <- 1-det(cor(data_num))
coef_dep_selector1 <- 1-det(cor(data_num_selector1))
coef_dep_selector2 <- 1-det(cor(data_num_selector2))
El coeficiente de dependencia del total de los datos es 0.9941683, mientras que para los datos etiquetados como Selector == 1 es 0.9917558 y para los datos etiquetados con Selector == 2 es 0.9998612.
Realizar un análisis de componentes principales y estudiar la calidad de la representación por diversos criterios según el número de dimensiones.
Primeramente, debemos tener en cuenta que el análisis de componentes principales (en adelante PCA) puede hacerse a partir de la matriz de covarianzas o la matriz de correlaciones.
Normalmente, si las distintas variables están medidas en escalas muy diferente y la varianza de las variables difiere mucho una de otra, no es recomendable que se ejecute el PCA a partir de la matriz de covarianzas, pues las primeras componentes principales van a estar definidas en base a una o pocas variables con una gran varianza.
A continuación observamos las varianzas de cada variable (presentes en la diagonal de la matriz de covarianzas):
diag(var(data_num)) %>% round(3) %>% knitr::kable(caption = "Tabla 1.3. Varianzas de las distintas variables de nuestro set de datos.")
| x | |
|---|---|
| Age | 263.146 |
| TB | 38.784 |
| DB | 7.933 |
| Alkphos | 59322.381 |
| Sgpt | 33555.955 |
| Sgot | 84013.042 |
| TP | 1.176 |
| ALB | 0.631 |
| A.G | 0.102 |
Como se puede observar en la Tabla 1.3, las varianzas de las distintas variables son muy diferentes, siendo para algunas variables (i.e. Alkphos) mucho mayor que para otras (i.e. TP), por lo tanto, lo más adecuado en este caso sería usar la matriz de correlaciones para hacer el PCA.
Para buscar las componentes principales, usaremos la función princomp(), a la qual debemos dar nuestro set de datos (con únicamente las variables numéricas sin NA) e indicar que queremos que haga el cálculo en base a la matriz de correlaciones (cor=T).
pca <- princomp(x = data_num, cor = T)
pca_vars <- pca$sdev**2
pca_sum <- summary(pca)
pca_sum
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 1.6594914 1.4236792 1.1685377 0.9787747 0.91905411
## Proportion of Variance 0.3059902 0.2252069 0.1517200 0.1064444 0.09385116
## Cumulative Proportion 0.3059902 0.5311971 0.6829171 0.7893616 0.88321272
## Comp.6 Comp.7 Comp.8 Comp.9
## Standard deviation 0.81632677 0.45105044 0.35470439 0.235445205
## Proportion of Variance 0.07404327 0.02260517 0.01397947 0.006159383
## Cumulative Proportion 0.95725598 0.97986115 0.99384062 1.000000000
prop_var_pca <- data.frame(pc = c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9"),
prop_var = c(0.3059902, 0.2252069, 0.1517200, 0.1064444, 0.09385116,
0.07404327, 0.02260517, 0.01397947, 0.006159383),
acum_var = c(0.3059902, 0.5311971, 0.6829171, 0.7893616, 0.88321272,
0.95725598, 0.97986115, 0.99384062, 1.000000000))
ggplot(prop_var_pca) +
geom_col(aes(x = pc, y = prop_var, group = 1, fill = pc), position = "dodge", show.legend = F) +
geom_point(aes(x = pc, y = acum_var/3,), show.legend = F) +
geom_line(aes(x = pc, y = acum_var/3, group = 1), show.legend = F) +
scale_y_continuous(name = expression("Explained variance (proportion)"),
sec.axis = sec_axis(~ . * 3 , name = "Cumulative variance (proportion)"), limits = c(0, 0.35)) +
xlab("Principal component") +
ggtitle("Variance explained by each principal component")
Figura 1.3. Proporción de la varianza explicada por cada componente principal (barras, eje izquierdo) junto con la proporción acumulada (puntos y líneas, eje derecho). Este tipo de gráfico a veces se llama ‘scree plot’
¿Cuántas componentes principales debemos escoger? Para escoger con cuantas componentes principales nos quedamos para continuar el análisis, podemos usar varios criterios:
Criterio de porcentaje: este criterio se basa en definir un porcentaje mínimo, por ejemplo 80%, de la varianza que queremos explicar con nuestras componentes principales:
En el caso de definir un porcentaje mínimo del 80%, nos quedaremos con las primeras 5 componentes principales, ya que así pasamos del 80% de varianza explicada, concretamente tenemos un 88.3% (Figura 1.3).
Sin embargo, en algunos casos, el punto de corte basado en el porcentaje se define cuando, al considerar una componente principal más, no se aumenta mucho en el porcentaje de varianza explicada. En otras palabras, hay una estabilización en el porcentaje de varianza explicada. Así pues, fijándonos en la Figura 1.3, se puede ver como de la componente principal 5 a la componente principal 6 aún hay un aumento considerable de la varianza explicada, mientras que a partir de la 6, la proporción de varianza acumulada aumenta muy poco con cada componente que añadimos, por lo que tendríamos razones para coger hasta 6 componentes principales.
Criterio de Kaiser: al usar la matriz de correlaciones para obtener las componentes principales, estamos asumiendo que las variables observables tienen varianza 1, por lo que una componente principal con una varianza por encima de 1 explica más variación de los datos que cualquier variable observable, mientras que una componente principal con varianza menor que 1 explica una menos variación que una de las variables observables y, por lo tanto, la descartaríamos para el análisis.
Siguiendo este criterio, escogeríamos las 3 primeras componentes principales (desviación estándard de 1.6594914, 1.4236792 y 1.1685377 respectivamente), ya que son las únicas con una desviación estándard mayor que 1 y, por ende, varianza mayor que 1. En este caso, con 3 componentes principales, explicaríamos un 68% de la varianza de los datos.
En algunos casos, se considera que el punto de corte debe ser 0.7, por lo que en este caso, nos quedaríamos con 5 componentes principales, dado que la 6a componente principal tiene una varianza de 0.6663894.
Calidad de la representación:
Para ver la calidad de la respresentación podemos usar la función get_pca_var() del paquete factoextra, que nos devuelve los resultados (coordenadas, cosenos cuadrados y contribuciones) de las variables en el PCA. En este caso nos interesan los cosenos cuadrados (cos2), que representan la calidad de la representación de las variables y es calculado como las coordenadas (coords) al cuadrado. En este caso, un valor alto de cos2 nos indica una buena representación de la variable en la dimensión o componente principal en cuestión.
pca_var <- factoextra::get_pca_var(res.pca = pca)
corrplot::corrplot(pca_var$cos2, is.corr=FALSE, method = "color", addgrid.col = "gray50", addCoef.col = "black",
title = "Cos2 of variables in each dimension or PC.")
Figura 1.4. Cos2 (calidad de la representación) de las variables en cada dimensión o componente principal. Un valor alto de ‘cos2’ indica que la variable tiene una buena representación en la dimensión en cuestión.
En la Figura 1.4 podemos ver la calidad de la representación de cada variable en cada dimensión: a un valor más alto de cos2 mejor representación tiene dicha variable. Por ejemplo, en la primera dimensión (o componente principal), las variables mejor representadas son la TB, la DB, la ALB y en menor medida la A.G. Hay algunas variables que tienen la mayor parte de su representación en las dos o tres primeras componentes principales, como la ALB o la Spgt, mientras que otras como Age se encuentran en la quinta dimensión y otras repartidas en las primeras componentes y las últimas (A.G).
A parte del mapa de calor de la Figura 1.4, también podemos visualizar un gráfico de barras donde se nos da la calidad de representación de las variables en un número determinado de dimensiones. Por ejemplo, en la Figura 1.5, vemos cuatro gráficos de barras indicandonos la calidad de la representación total de las variables en las 2, 3, 5 o 6 primeras componentes principales.
themes <- theme(axis.title.y = element_text(size = 9))
cos1_2 <- fviz_cos2(pca, choice = "var", axes = 1:2, fill = "darkblue", color = "darkblue") + themes
cos1_3 <- fviz_cos2(pca, choice = "var", axes = 1:3, fill = "blue2", color = "blue2") + themes
cos1_5 <- fviz_cos2(pca, choice = "var", axes = 1:5, fill = "cadetblue", color = "cadetblue") + themes
cos1_6 <- fviz_cos2(pca, choice = "var", axes = 1:6, fill = "lightblue", color = "lightblue") + themes
(cos1_2+cos1_3)/(cos1_5+cos1_6)
Figura 1.5. Calidad de represntacion de las variables en las 2, 3, 5 o 6 primeras componentes principales.
Algo que se ve claramente es que con las dos primeras dimensiones, hay una gran diferencia de calidad de representación entre la ALB (la mejor representada) y las otras variables, que mayoritariamanete se encuentran por encima de 0.5, salvo Age y Alkphos.
Si aumentamos el número de dimensiones y cogemos las tres primeras, lo que hemos considerado con el criterio del Kaiser (con un corte de 1), la calidad de la representación de la mayoría de las variables aumenta (salvo por A.G, Age y Alkphos).
Si aumentamos las dimensiones a las 5 primeras (que explican casi el 90% de la variablidad y hemos considerado con el criterio de Kaiser con un corte de 0.7), tenemos que casi todas las variables (salvo A.G) estan en un valor de cos2 mayor que 0.75 indicando una muy buena representación.
Finalmente, con las seis primeras componentes (que hemos escogido según el criterio gráfico ), todas las variables tienen una calidad de representación muy alta, con un cos2 cerca de 1.
Alternativamente, también se pueden visualizar los datos mediante las coordenadas por pares de dimensiones con la función fviz_pca_var(). Esta función no nos da exactamente la calidad de la representación, sinó que representa las variables usando las coordenadas para cada dimensión y visualizando las dimensiones de 2 en 2, de esta manera podemos observar las correlaciones entre variables en dichas dimensiones. Además, si consideramos un gradiente de colores en base al valor de cos2, podemos observar la calidad de la representación en las dimensiones observadas, tal y como se ha hecho en la Figura S1.2. Así mismo, este tipo de gráficos o los gráficos de barras del estilo de la Figura 1.5 pueden usarse para ver la contribución de cada variable a cada una de las dimensiones (funciones fviz_pca_contrib() y fviz_contrib()).
Cabe destacar que nos hemos centrado en la calidad de representación de las variables, pero de la misma forma, podríamos estudiar la calidad de la representación de los individuos si usaramos la funciones análogas a las usadas (i.e. fviz_cos2(pca, choice = "ind", ...), etc).
Interpretar los dos primeros ejes principales mediante el gráfico de correlaciones con las variables originales.
El gráfico de correlaciones se obtiene con la función fviz_pca_var(), como se ha hecho en las figuras S1.2 y 1.6. Con este gráfico podemos ver las correlaciones entre las distintas variables en las dimensiones indicadas (en este caso PC1 y PC2):
fviz_pca_var(pca, col.var = "contrib", gradient.cols =c("#00AFBB", "#E7B800", "#FC4E07"))
Figura 1.6. Gráfico de correlaciones de las variables en las dos primeras componentes principales. El gradiente de color indica la contribución de cada variable a las dimensiones correspondientes.
En la Figura 1.6 vemos como hay algunas variables muy correlacionadas en las primeras dos dimensiones. Por ejemplo, TB y DB tienen una correlación casi perfecta en las componenentes principales 1 y 2, así como Sgot y Sgpt. En cambio, vemos coo la variable ALB está inversamente correalcionada con la variable Age, aunque vasándonos en su contribución y su calidad de representación, no podemos hacer predicciones sobre la variable Age en las componentes 1 y 2.
Representar los pacientes en un gráfico de dos dimensiones con puntos distintos según el grupo y el género.
pca <- princomp(data_filt %>% dplyr::select(everything(), -Gender, -Selector), cor = T)
# mutate datafilt to create a new variable mixing gender and selector
group <- data_filt %>% mutate(SelGen = paste(Gender, "-", Selector))
fviz_pca_ind(pca,
geom.ind = "point",
fill.ind = group$SelGen,
pointshape = 21, pointsize = 2,
legend.title = list(fill = "Gender", color = "Selector"))
Figura 1.7. Gráfico en als dimensiones 1 (eje X) y 2 (eje Y) de los pacientes segun el grupo (Selector) y género (Gender).
¿Alguien se atreve con el gráfico en tres dimensiones?
El gráfico de el PCA en tres dimensiones se puede llevar a cabo mediante el paquete pca3d.
library(pca3d)
pc <- princomp(data_num, cor=TRUE, scores=TRUE)
gr <- data_filt$Gender
pca3d(pca = pc, group = gr, components = c(1,2,3), show.scale = T, show.plane = T, legend = "topleft", title = "3D-PCA plot on the first 3 PCs. Gender coloring.")
snapshotPCA3d(file="figs/PCA3d_gender.png")
gr <- data_filt$Selector
pca3d(pca = pc, group = gr, components = c(1,2,3), show.scale = T, show.plane = T, legend = "topleft", title = "3D-PCA plot on the first 3 PCs. Selector coloring.")
snapshotPCA3d(file="figs/PCA3d_selector.png")
gr <- data_filt %>% mutate(SelGen = paste(Gender, "-", Selector))
gr <- gr$SelGen
pca3d(pca = pc, group = gr, components = c(1,2,3), show.scale = T, show.plane = T, legend = "topleft", title = "3D-PCA plot on the first 3 PCs. Selector coloring.")
snapshotPCA3d(file="figs/PCA3d_SelGen.png")
knitr::include_graphics("figs/PCA3d_SelGen.png")
Figura 1.8. Gráfico 3D de los pacientes en las tres primeras componentes. Los puntos están coloresados según el género y el grupo de cada uno.
En la Figura S1.3 se puede observar un gráfico similar, pero usando la función scatterplot3d() del paquete con el mismo nombre.
Una de las finalidades de disponer de marcadores para la enfermedad es determinar si ciertas variables se encuentran alteradas en los pacientes afectados. Para ello vamos a utilizar los datos del ejercicio 1, para investigar si hay diferencias entre los dos grupos.
a. Realizar un boxplot múltiple que permita comparar visualmente todas las variables numéricas, separadas por controles y pacientes.
Para hacer un boxplot múltiple, primero de todo usaremos la dunción melt() del paquete reshape2 y después usaremos el paquete ggplot2 (funciones ggplot() y geom_boxplot()). Finalmente usaremos la función facet_wrap() para separar los valores de la variable Selector en dos paneles o, por otro lado, también podremos separar todas las variables en distintos paneles para hacer más fácil la comparación entre pacientes y controles por cada variable.
data_filt %>%
dplyr::mutate(Selector = Selector %>% recode("1" = "Pacients","2" = "Controls")) %>%
melt() %>%
ggplot(aes(x=Selector, y=value, fill=Selector)) + geom_boxplot(outlier.colour = "Gray50") +
facet_wrap(~variable, scales = "free_y") +
stat_compare_means(method = "wilcox.test", paired = F,
label = "p.format", label.y.npc = 0.85, label.x.npc = 0.3) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 0.95),
legend.position = "none") +
ggtitle("Variables of our study separated in patients and control subjects.") +
ylab("") + xlab("")
Figura 2.1. Gráficos de cajas de las variables numéricas presentes en los datos del ejercicio 1, separados según su condición (‘Selector’) de pacientes o controles. El p-valor mostrado ha sido obtenido mediante un test de wilcoxon (no pareado) e indica si las medias de pacientes y controles para la variable correspondiente son significativamente distintas o no.
La Figura 2.1 nos permite comparar fácilmente muchas de las variables entre los pacientes (Selector = 1) y los controles (Selector = 2) de nuestro estudio:
La edad de los pacientes (mediana \(\simeq\) 46) es ligeramente mayor que los controles (mediana \(\simeq\) 41).
Aunque es dificil de ver a simple vista la DB y la TB de los pacientes (medianas: DB \(\simeq\) 0.5, TB \(\simeq\) 1.4) es ligeramente superior a la de los controles (medianas: DB \(\simeq\) 0.2, TB \(\simeq\) 0.8).
Por lo que respeta a las variables TP, ALB y A.G, se ve una ligera diferencia, siendo sus valores mayores en el caso del grupo control (medianas: TP \(\simeq\) 6.6, ALB \(\simeq\) 3.4, A.G \(\simeq\) 1), que en el grupo de los pacientes (medianas: TP \(\simeq\) 6.55, ALB \(\simeq\) 3, A.G \(\simeq\) 0.9). Sin embargo, esta diferencia no es significativa para la variable TP.
Finalmente están las variables Alkphos, Sgpt y Sgot, las medias són significativamente diferentes y con unos valores mayores en el caso de los pacientes (medianas: Alkphos \(\simeq\) 229, Sgpt \(\simeq\) 51, Sgot \(\simeq\) 53) que de los sujetos del grupo control (medianas: Alkphos \(\simeq\) 187, Sgpt \(\simeq\) 28, Sgot \(\simeq\) 29).
En la Figura S2.1 se puede ver un gráfico con los mismos datos separados en dos paneles, uno para pacientes y otro para controles.
Si suponemos que se verifican todas las condiciones, realizar un MANOVA de un factor, utilizando la variable Selector, para decidir si hay diferencias en las variables marcadores
Antes de realizar el test MANOVA debemos hacer varias suposciciones:
mshapiro_test().mshapiro_test(data = data_num)
## # A tibble: 1 x 2
## statistic p.value
## <dbl> <dbl>
## 1 0.124 2.55e-45
Homogeneidad de las varianzas: esto se puede llevar a cabo con un test de Levene (varianzas) o con un test de Box (covarianzas). Más adelante veremos que las varianzas y covarianzas entre grupos no son homogeneas.
…
Como hemos visto, los datos no cumplen una distribución normal, ya que el p-valor obtenido del test de Shapiro-Wilk, es mucho menor que el nivel de significancia de 0.05, con lo que no podemos asumir la hipòtesis nula de la distribución normal. Así pues, tendríamos que usar un test no paramétrico homólogo al MANOVA.
Sin embargo, asumiendo que todas las condiciones mencionadas se cumplen, para hacer el test MANOVA en un factor (Selector) podemos usar la función manova() que admite como argumentos los mismos que la función aov para el ANOVA.
attach(data_filt)
manova <- manova(cbind(A.G, ALB, TP, Sgot, Sgpt, Alkphos, DB, TB, Age) ~ Selector, data = data_filt)
detach()
summary(manova)
## Df Pillai approx F num Df den Df Pr(>F)
## Selector 1 0.11654 8.34 9 569 9.984e-12 ***
## Residuals 577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Segun el test MANOVA que hemos realizado, sí que hay diferencias significativas en las variables marcadores.
summary.aov(manova)
## Response A.G :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 1.571 1.57107 15.775 8.037e-05 ***
## Residuals 577 57.465 0.09959
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response ALB :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 9.31 9.3118 15.114 0.0001129 ***
## Residuals 577 355.48 0.6161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response TP :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 0.77 0.76831 0.6527 0.4195
## Residuals 577 679.22 1.17715
##
## Response Sgot :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 1119477 1119477 13.616 0.0002455 ***
## Residuals 577 47440061 82218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sgpt :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 516055 516055 15.772 8.049e-05 ***
## Residuals 577 18879287 32720
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Alkphos :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 1152846 1152846 20.075 8.982e-06 ***
## Residuals 577 33135491 57427
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response DB :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 278.1 278.087 37.255 1.903e-09 ***
## Residuals 577 4307.0 7.464
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response TB :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 1087.2 1087.15 29.408 8.633e-08 ***
## Residuals 577 21330.3 36.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Age :
## Df Sum Sq Mean Sq F value Pr(>F)
## Selector 1 2697 2697.09 10.416 0.00132 **
## Residuals 577 149401 258.93
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Concretamente, con un nivel de significancia de 0.05, todas las variables numéricas salvo TP muestran diferencias significativas entre los dos grupos.
Como además tenemos el factor Gender podemos considerar un modelo general con todos los factores (sin interacciones) y el modelo más simple con solo el factor principal y contrastar con la función anova() si el modelo simple explica la variación de los datos.
Para hacer el modelo general con los dos factores (Selector y Gender) usaremos la función manova(), como hicimos en el apartado B, donde generamos el modelo simple con un solo factor (Selector). Luego para comparar los dos modelos, podemos usar la función summary() o la función anova() y ver que factores explican significativamente la variación de los datos.
attach(data_filt)
manova_general <- manova(cbind(A.G, ALB, TP, Sgot, Sgpt, Alkphos, DB, TB, Age) ~ Selector+Gender, data = data_filt)
detach()
Con el model simple de un factor:
anova(manova)
## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98774 5095.5 9 569 < 2.2e-16 ***
## Selector 1 0.11654 8.3 9 569 9.984e-12 ***
## Residuals 577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vemos como la variable Selector explica significativamente la variación de los datos.
Con el model general con los dos factores:
anova(manova_general)
## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98775 5086.9 9 568 < 2.2e-16 ***
## Selector 1 0.11717 8.4 9 568 8.788e-12 ***
## Gender 1 0.02898 1.9 9 568 0.05186 .
## Residuals 576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vemos como la variable Selector explica significativamente la variación de los datos, mientras que la variable Gender no es significativa (usando un nivel de significancia de 0.05, pero sí si usamos de 0.1).
Así pues, como en el modelo general el único factor que explica la variación de los datos es Selector, podemos concluir que con el modelo simple tenemos suficiente para explicar dicha variación.
Comparar las matrices de covarianzas de los dos grupos del factor Selector con el test de la razón de verosimilitudes y con el test M de Box
Las matrices de covarianzas se pueden encontrar en la Figura S2.1 para Selector = 1 y en la Figura S2.2 para Selector = 2.
Test de razón de verosimilitudes:
Test M de Box: el test M de Box, o prueba de box permite comparar las matrices de covarianzas entre dos grupos para ver si estas son iguales. Es un test muy sensible y debido a esto, se recomienda un límite de significancia bastante bajo (i.e. 0.001). También hay que tener en cuenta, que en caso de falta de normalidad de los datos, el resultado del test de Box puede ser significativo debido a esta falta de normalidad y no a la diferencia entre las matrices de covarianzas. Para llevar a cabo el test, se puede usar la función boxM() del paquete heplots.
boxM_res <- heplots::boxM(data_num, data_filt$Selector)
summary(boxM_res)
## Summary for Box's M-test of Equality of Covariance Matrices
##
## Chi-Sq: 2496.705
## df: 45
## p-value: < 2.2e-16
##
## log of Covariance determinants:
## 1 2 pooled
## 38.00520 16.17060 36.22022
##
## Eigenvalues:
## 1 2 pooled
## 1 1.488917e+05 2.023218e+04 1.074514e+05
## 2 6.973537e+04 1.505127e+03 5.519672e+04
## 3 1.346050e+04 2.889142e+02 9.726111e+03
## 4 2.378558e+02 2.647388e+02 2.532873e+02
## 5 5.525106e+01 1.527346e+00 3.979737e+01
## 6 2.193000e+00 7.650083e-01 1.573901e+00
## 7 1.542026e+00 1.544093e-01 1.548583e+00
## 8 2.083380e-01 5.788200e-03 1.973065e-01
## 9 2.474629e-02 4.332887e-03 1.921571e-02
##
## Statistics based on eigenvalues:
## 1 2 pooled
## product 3.202195e+16 1.053903e+07 5.373311e+15
## sum 2.323846e+05 2.229342e+04 1.726707e+05
## precision 2.158108e-02 2.427158e-03 1.711762e-02
## max 1.488917e+05 2.023218e+04 1.074514e+05
Como se puede observar, con la prueba de Box, obtenemos que las matrices de covarianzas de los dos grupos son significativamente distintas. Aunque bien podría deberse a la falta de normalidad de los datos.
Otra opción para el apartado anterior habría sido comparar la variabilidad con un test de Levene múltiple. ¿Este test contrasta lo mismo que los anteriores?
El test de Levene es un test diseñado para el análisis de varianzas. Está diseñado para análisis univariante de la varianza (ANOVA), aunque se puede extrapolar a multivariante, pero con este contrastamos la homogeneidad de cada variable (igualdad de varianzas) entre los grupos, mientras con los dos tests efectuados en el apartado D, miramos la varianza y las covarianzas de todas las variables entre los grupos estudiados (Selector = 1, Selector = 2).
library(car)
attach(data_filt)
levene_test(data = data_filt, formula = A.G+ALB+TP+Sgot+Sgpt+Alkphos+DB+TB+Age ~ Selector)
## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 577 23.5 0.00000158
detach()
Los experimentadores dudan de la normalidad de las variables observadas. Estudiar el caso.
Gràficamente, podemos usar la función mvn() para hacer los gráficos univariantes (i.e. histogramas) y multivariante (i.e. qqplot):
mvn_res <- MVN::mvn(data = data_num_selector1, univariatePlot = "histogram")
Figura 2.2a. Histogramas de para cada variable en los datos de pacientes (Selector = 1).
mvn_res <- MVN::mvn(data = data_num_selector2, univariatePlot = "histogram")
Figura 2.2b. Histogramas de para cada variable en los datos de sujetos control (Selector = 2).
mvn_res <- MVN::mvn(data = data_num_selector1, multivariatePlot = "qq")
Figura 2.3a. QQplot de las distancias de Mahalanobis en los datos de pacientes (Selector = 1)
mvn_res <- MVN::mvn(data = data_num_selector2, multivariatePlot = "qq")
Figura 2.3b. QQplot de las distancias de Mahalanobis en los datos de sujetos control (Selector = 2)
A nivel gráfico, observando las Figuras 2.2a y 2.2b se podría intuir una normalidad univariante para las variables Age y TP en ambos grupos, mientras que ALB parece tener una distribución similar a la normal en el sólo grupo de los pacientes y las otras variable no se parecen para nada a una distribución normal, pues són altamente asimétricas en ambos grupos.
En cuanto a los gráficos (qqplots) multivariantes (Figuras 2.3a y 2.3b), estos nos dicen que los datos no siguen una normalidad multivariante.
Sin embargo, no podemos fiarnos solo de los métodos gráficos y se tienen que usar otros test estadísticos, como Shapiro-Wilk, que nos servirá tanto a nivel univariante como multivariante. * Shapiro-Wilk univariante:
shapiro_selector1 <- shapiro_test(data_num_selector1, vars = colnames(data_num_selector1))
shapiro_selector2 <- shapiro_test(data_num_selector2, vars = colnames(data_num_selector1))
paste("Resultado del test Shapiro-Wilk para el grupo Selector 1:")
## [1] "Resultado del test Shapiro-Wilk para el grupo Selector 1:"
shapiro_selector1
## # A tibble: 9 x 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 A.G 0.928 3.31e-13
## 2 Age 0.990 6.28e- 3
## 3 ALB 0.993 4.59e- 2
## 4 Alkphos 0.624 4.29e-29
## 5 DB 0.607 1.39e-29
## 6 Sgot 0.320 2.16e-36
## 7 Sgpt 0.374 2.61e-35
## 8 TB 0.529 9.54e-32
## 9 TP 0.991 1.32e- 2
paste("Como se puede ver, todas las variables muestran un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 1.")
## [1] "Como se puede ver, todas las variables muestran un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 1."
paste("")
## [1] ""
paste("Resultado del test Shapiro-Wilk para el grupo Selector 2:")
## [1] "Resultado del test Shapiro-Wilk para el grupo Selector 2:"
shapiro_selector2
## # A tibble: 9 x 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 A.G 0.969 9.15e- 4
## 2 Age 0.983 4.33e- 2
## 3 ALB 0.980 1.72e- 2
## 4 Alkphos 0.495 1.42e-21
## 5 DB 0.506 2.33e-21
## 6 Sgot 0.635 1.28e-18
## 7 Sgpt 0.686 2.54e-17
## 8 TB 0.500 1.75e-21
## 9 TP 0.989 2.21e- 1
paste("Como se puede ver, todas las variables muestran menos TP un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 2.")
## [1] "Como se puede ver, todas las variables muestran menos TP un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 2."
mshapiro_selector1 <- mshapiro_test(data = data_num_selector1)
mshapiro_selector2 <- mshapiro_test(data = data_num_selector2)
paste("El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 1 es", mshapiro_selector1$p.value, "por lo que no podemos considerar que este grupo siga una distribución normal multivariante.")
## [1] "El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 1 es 2.00419808624749e-39 por lo que no podemos considerar que este grupo siga una distribución normal multivariante."
paste("El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 2 es", mshapiro_selector2$p.value, "por lo que no podemos considerar que este grupo siga una distribución normal multivariante.")
## [1] "El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 2 es 1.16805887608274e-21 por lo que no podemos considerar que este grupo siga una distribución normal multivariante."
mvn() del paquete MVN. La misma función permite el estudio univariante y también hace los gráficos.mvn_res <- MVN::mvn(data = data_num)
mvn_res$multivariateNormality %>% knitr::kable(caption = "Tabla 2.1. Resultado de la función 'mvn()' para la normalidad multivariante.")
| Test | Statistic | p value | Result |
|---|---|---|---|
| Mardia Skewness | 71910.4739009141 | 0 | NO |
| Mardia Kurtosis | 788.563271538733 | 0 | NO |
| MVN | NA | NA | NO |
mvn_res$univariateNormality %>% knitr::kable(caption = "Tabla 2.2. Resultado de la función 'mvn()' para la normalidad univariante.")
| Test | Variable | Statistic | p value | Normality |
|---|---|---|---|---|
| Shapiro-Wilk | Age | 0.9920 | 0.0033 | NO |
| Shapiro-Wilk | TB | 0.4613 | <0.001 | NO |
| Shapiro-Wilk | DB | 0.5312 | <0.001 | NO |
| Shapiro-Wilk | Alkphos | 0.5859 | <0.001 | NO |
| Shapiro-Wilk | Sgpt | 0.3286 | <0.001 | NO |
| Shapiro-Wilk | Sgot | 0.2811 | <0.001 | NO |
| Shapiro-Wilk | TP | 0.9918 | 0.0029 | NO |
| Shapiro-Wilk | ALB | 0.9925 | 0.0053 | NO |
| Shapiro-Wilk | A.G | 0.9465 | <0.001 | NO |
d2 <- Moutlier(data_num, quantile = 0.975, plot = F)
n <- nrow(data_num)
p <- ncol(data_num)
quantiles <- qchisq((1:n)/(n+1),p)
plot(quantiles, sort(d2$md), ylab = "Ordered classical Mahalanobis distances",xlab = "Chi-square quantile")
Figura 2.4. QQplot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis clàsicas.
plot(quantiles, sort(d2$rd), ylab = "Ordered robust Mahalanobis distances", xlab = "Chi-square quantile")
Figura 2.5. QQplot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis robustas.
Según cual haya sido el resultado de los apartados anteriores y dado que son únicamente dos grupos, tal vez se deba utilizar un test de permutaciones con la T2 de Hotelling (Estudiar la función hotelling.test() del paquete Hotelling)
cor_list_1d %>% knitr::kable(align = "c", caption = "Tabla S1.1. Lista de correlaciones dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos.")
| Variable 1 | Variable 2 | Pearson’s correlation | Absolute Pearson’s |
|---|---|---|---|
| Age | Age | 1.0000000 | 1.0000000 |
| TB | TB | 1.0000000 | 1.0000000 |
| DB | DB | 1.0000000 | 1.0000000 |
| Alkphos | Alkphos | 1.0000000 | 1.0000000 |
| Sgpt | Sgpt | 1.0000000 | 1.0000000 |
| Sgot | Sgot | 1.0000000 | 1.0000000 |
| TP | TP | 1.0000000 | 1.0000000 |
| ALB | ALB | 1.0000000 | 1.0000000 |
| A.G | A.G | 1.0000000 | 1.0000000 |
| DB | TB | 0.8744810 | 0.8744810 |
| TB | DB | 0.8744810 | 0.8744810 |
| Sgot | Sgpt | 0.7918621 | 0.7918621 |
| Sgpt | Sgot | 0.7918621 | 0.7918621 |
| ALB | TP | 0.7831122 | 0.7831122 |
| TP | ALB | 0.7831122 | 0.7831122 |
| A.G | ALB | 0.6896323 | 0.6896323 |
| ALB | A.G | 0.6896323 | 0.6896323 |
| ALB | Age | -0.2642109 | 0.2642109 |
| Age | ALB | -0.2642109 | 0.2642109 |
| Sgot | DB | 0.2570224 | 0.2570224 |
| DB | Sgot | 0.2570224 | 0.2570224 |
| Sgot | TB | 0.2373231 | 0.2373231 |
| TB | Sgot | 0.2373231 | 0.2373231 |
| A.G | TP | 0.2348872 | 0.2348872 |
| TP | A.G | 0.2348872 | 0.2348872 |
| A.G | Alkphos | -0.2341665 | 0.2341665 |
| Alkphos | A.G | -0.2341665 | 0.2341665 |
| Alkphos | DB | 0.2340076 | 0.2340076 |
| DB | Alkphos | 0.2340076 | 0.2340076 |
| Sgpt | DB | 0.2331801 | 0.2331801 |
| DB | Sgpt | 0.2331801 | 0.2331801 |
| ALB | DB | -0.2284092 | 0.2284092 |
| DB | ALB | -0.2284092 | 0.2284092 |
| ALB | TB | -0.2220866 | 0.2220866 |
| TB | ALB | -0.2220866 | 0.2220866 |
| A.G | Age | -0.2164083 | 0.2164083 |
| Age | A.G | -0.2164083 | 0.2164083 |
| Sgpt | TB | 0.2133755 | 0.2133755 |
| TB | Sgpt | 0.2133755 | 0.2133755 |
| A.G | TB | -0.2062672 | 0.2062672 |
| TB | A.G | -0.2062672 | 0.2062672 |
| Alkphos | TB | 0.2057392 | 0.2057392 |
| TB | Alkphos | 0.2057392 | 0.2057392 |
| A.G | DB | -0.2001247 | 0.2001247 |
| DB | A.G | -0.2001247 | 0.2001247 |
| TP | Age | -0.1862481 | 0.1862481 |
| Age | TP | -0.1862481 | 0.1862481 |
| Sgot | Alkphos | 0.1665800 | 0.1665800 |
| Alkphos | Sgot | 0.1665800 | 0.1665800 |
| ALB | Alkphos | -0.1634186 | 0.1634186 |
| Alkphos | ALB | -0.1634186 | 0.1634186 |
| Sgpt | Alkphos | 0.1247767 | 0.1247767 |
| Alkphos | Sgpt | 0.1247767 | 0.1247767 |
| Sgpt | Age | -0.0877992 | 0.0877992 |
| Age | Sgpt | -0.0877992 | 0.0877992 |
| ALB | Sgot | -0.0849146 | 0.0849146 |
| Sgot | ALB | -0.0849146 | 0.0849146 |
| Alkphos | Age | 0.0788784 | 0.0788784 |
| Age | Alkphos | 0.0788784 | 0.0788784 |
| A.G | Sgot | -0.0700398 | 0.0700398 |
| Sgot | A.G | -0.0700398 | 0.0700398 |
| TP | Sgpt | -0.0424321 | 0.0424321 |
| Sgpt | TP | -0.0424321 | 0.0424321 |
| ALB | Sgpt | -0.0286575 | 0.0286575 |
| Sgpt | ALB | -0.0286575 | 0.0286575 |
| TP | Alkphos | -0.0270620 | 0.0270620 |
| Alkphos | TP | -0.0270620 | 0.0270620 |
| TP | Sgot | -0.0257510 | 0.0257510 |
| Sgot | TP | -0.0257510 | 0.0257510 |
| Sgot | Age | -0.0204989 | 0.0204989 |
| Age | Sgot | -0.0204989 | 0.0204989 |
| TB | Age | 0.0110004 | 0.0110004 |
| Age | TB | 0.0110004 | 0.0110004 |
| TP | TB | -0.0079059 | 0.0079059 |
| TB | TP | -0.0079059 | 0.0079059 |
| DB | Age | 0.0067843 | 0.0067843 |
| Age | DB | 0.0067843 | 0.0067843 |
| A.G | Sgpt | -0.0023750 | 0.0023750 |
| Sgpt | A.G | -0.0023750 | 0.0023750 |
| TP | DB | 0.0000327 | 0.0000327 |
| DB | TP | 0.0000327 | 0.0000327 |
par_cor_1e_list <- par_cor_1e$estimate %>%
reshape2::melt() %>%
dplyr::mutate(abs_pearson = sqrt(value^2)) %>%
dplyr::arrange(desc(abs_pearson)) %>%
magrittr::set_colnames(c("Variable 1", "Variable 2", "Pearson's correlation", "Absolute Pearson's"))
par_cor_1e_list %>% knitr::kable(align = "c", caption = "Tabla 1.2. Lista de correlaciones parciales dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos.")
| Variable 1 | Variable 2 | Pearson’s correlation | Absolute Pearson’s |
|---|---|---|---|
| Age | Age | 1.0000000 | 1.0000000 |
| TB | TB | 1.0000000 | 1.0000000 |
| DB | DB | 1.0000000 | 1.0000000 |
| Alkphos | Alkphos | 1.0000000 | 1.0000000 |
| Sgpt | Sgpt | 1.0000000 | 1.0000000 |
| Sgot | Sgot | 1.0000000 | 1.0000000 |
| TP | TP | 1.0000000 | 1.0000000 |
| ALB | ALB | 1.0000000 | 1.0000000 |
| A.G | A.G | 1.0000000 | 1.0000000 |
| ALB | TP | 0.8971925 | 0.8971925 |
| TP | ALB | 0.8971925 | 0.8971925 |
| DB | TB | 0.8401910 | 0.8401910 |
| TB | DB | 0.8401910 | 0.8401910 |
| A.G | ALB | 0.8246399 | 0.8246399 |
| ALB | A.G | 0.8246399 | 0.8246399 |
| Sgot | Sgpt | 0.7829975 | 0.7829975 |
| Sgpt | Sgot | 0.7829975 | 0.7829975 |
| A.G | TP | -0.6865110 | 0.6865110 |
| TP | A.G | -0.6865110 | 0.6865110 |
| DB | ALB | -0.1910744 | 0.1910744 |
| ALB | DB | -0.1910744 | 0.1910744 |
| TP | DB | 0.1823035 | 0.1823035 |
| DB | TP | 0.1823035 | 0.1823035 |
| TP | Sgpt | -0.1656768 | 0.1656768 |
| Sgpt | TP | -0.1656768 | 0.1656768 |
| ALB | Sgpt | 0.1469404 | 0.1469404 |
| Sgpt | ALB | 0.1469404 | 0.1469404 |
| A.G | DB | 0.1412185 | 0.1412185 |
| DB | A.G | 0.1412185 | 0.1412185 |
| TP | Sgot | 0.1234029 | 0.1234029 |
| Sgot | TP | 0.1234029 | 0.1234029 |
| ALB | Sgot | -0.1171630 | 0.1171630 |
| Sgot | ALB | -0.1171630 | 0.1171630 |
| Alkphos | A.G | -0.1165394 | 0.1165394 |
| A.G | Alkphos | -0.1165394 | 0.1165394 |
| Alkphos | DB | 0.0971593 | 0.0971593 |
| DB | Alkphos | 0.0971593 | 0.0971593 |
| Age | Sgpt | -0.0957684 | 0.0957684 |
| Sgpt | Age | -0.0957684 | 0.0957684 |
| ALB | Age | -0.0820090 | 0.0820090 |
| Age | ALB | -0.0820090 | 0.0820090 |
| A.G | Sgpt | -0.0715968 | 0.0715968 |
| Sgpt | A.G | -0.0715968 | 0.0715968 |
| Alkphos | Sgot | 0.0700443 | 0.0700443 |
| Sgot | Alkphos | 0.0700443 | 0.0700443 |
| A.G | Sgot | 0.0591914 | 0.0591914 |
| Sgot | A.G | 0.0591914 | 0.0591914 |
| Sgpt | DB | 0.0579524 | 0.0579524 |
| DB | Sgpt | 0.0579524 | 0.0579524 |
| Sgot | Age | 0.0523929 | 0.0523929 |
| Age | Sgot | 0.0523929 | 0.0523929 |
| Age | Alkphos | 0.0461958 | 0.0461958 |
| Alkphos | Age | 0.0461958 | 0.0461958 |
| A.G | TB | -0.0293029 | 0.0293029 |
| TB | A.G | -0.0293029 | 0.0293029 |
| Age | DB | -0.0248805 | 0.0248805 |
| DB | Age | -0.0248805 | 0.0248805 |
| A.G | Age | -0.0228778 | 0.0228778 |
| Age | A.G | -0.0228778 | 0.0228778 |
| Alkphos | ALB | 0.0146719 | 0.0146719 |
| ALB | Alkphos | 0.0146719 | 0.0146719 |
| Alkphos | TB | -0.0135221 | 0.0135221 |
| TB | Alkphos | -0.0135221 | 0.0135221 |
| Sgot | TB | 0.0118750 | 0.0118750 |
| TB | Sgot | 0.0118750 | 0.0118750 |
| DB | Sgot | 0.0105666 | 0.0105666 |
| Sgot | DB | 0.0105666 | 0.0105666 |
| ALB | TB | -0.0087443 | 0.0087443 |
| TB | ALB | -0.0087443 | 0.0087443 |
| TB | TP | 0.0079615 | 0.0079615 |
| TP | TB | 0.0079615 | 0.0079615 |
| TP | Age | 0.0075072 | 0.0075072 |
| Age | TP | 0.0075072 | 0.0075072 |
| TB | Sgpt | 0.0061201 | 0.0061201 |
| Sgpt | TB | 0.0061201 | 0.0061201 |
| TB | Age | -0.0029639 | 0.0029639 |
| Age | TB | -0.0029639 | 0.0029639 |
| Sgpt | Alkphos | 0.0010541 | 0.0010541 |
| Alkphos | Sgpt | 0.0010541 | 0.0010541 |
| Alkphos | TP | -0.0002238 | 0.0002238 |
| TP | Alkphos | -0.0002238 | 0.0002238 |
cov_mat_selector1 %>% knitr::kable(align = "c", caption = "Tabla S2.1. Matriz de covarianzas del grupo de pacientes ('Selector' = 1).")
| Age | TB | DB | Alkphos | Sgpt | Sgot | TP | ALB | A.G | |
|---|---|---|---|---|---|---|---|---|---|
| Age | 246.187178 | -2.6722181 | -1.7319577 | 376.866407 | -436.8823385 | -243.458680 | -3.1780187 | -3.4973541 | -1.3308976 |
| TB | -2.672218 | 51.2423524 | 19.9290225 | 318.895747 | 279.9389462 | 509.423992 | 0.0609333 | -1.2364565 | -0.4557098 |
| DB | -1.731958 | 19.9290225 | 10.3203765 | 164.754420 | 137.6199711 | 247.996179 | 0.0670728 | -0.5629437 | -0.1916362 |
| Alkphos | 376.866407 | 318.8957469 | 164.7544198 | 72277.363091 | 5390.7891006 | 13157.741868 | -5.9794926 | -29.9190722 | -18.3494859 |
| Sgpt | -436.882339 | 279.9389462 | 137.6199711 | 5390.789101 | 45461.4641834 | 56759.007053 | -10.6640541 | -0.8768935 | 1.9592638 |
| Sgot | -243.458680 | 509.4239920 | 247.9961785 | 13157.741868 | 56759.0070534 | 114336.114960 | -7.4678598 | -18.6288346 | -5.9487430 |
| TP | -3.178019 | 0.0609333 | 0.0670728 | -5.979493 | -10.6640541 | -7.467860 | 1.2040283 | 0.6576987 | 0.0720592 |
| ALB | -3.497354 | -1.2364565 | -0.5629437 | -29.919072 | -0.8768935 | -18.628835 | 0.6576987 | 0.6200131 | 0.1702830 |
| A.G | -1.330898 | -0.4557098 | -0.1916362 | -18.349486 | 1.9592638 | -5.948743 | 0.0720592 | 0.1702830 | 0.1063755 |
cov_mat_selector2 %>% knitr::kable(align = "c", caption = "Tabla S2.2. Matriz de covarianzas del grupo control ('Selector' = 2).")
| Age | TB | DB | Alkphos | Sgpt | Sgot | TP | ALB | A.G | |
|---|---|---|---|---|---|---|---|---|---|
| Age | 291.0133038 | 0.2049335 | 0.1732816 | -190.695676 | -46.7998891 | -61.6452328 | -3.2686807 | -2.2266075 | -0.2056375 |
| TB | 0.2049335 | 1.0193178 | 0.5197982 | 80.923976 | 8.5153104 | 14.2240798 | -0.1654361 | -0.1452531 | -0.0472918 |
| DB | 0.1732816 | 0.5197982 | 0.2724257 | 41.683603 | 4.3908647 | 7.3869401 | -0.0794290 | -0.0732705 | -0.0248308 |
| Alkphos | -190.6956763 | 80.9239763 | 41.6836031 | 20030.119586 | 1341.8261641 | 1384.5409091 | -4.3997044 | -16.1198263 | -9.8257443 |
| Sgpt | -46.7998891 | 8.5153104 | 4.3908647 | 1341.826164 | 632.3328160 | 610.3452328 | 0.9814856 | 0.8766075 | 0.0663326 |
| Sgot | -61.6452328 | 14.2240798 | 7.3869401 | 1384.540909 | 610.3452328 | 1336.8645233 | -4.0711197 | -2.3125831 | 0.2007528 |
| TP | -3.2686807 | -0.1654361 | -0.0794290 | -4.399704 | 0.9814856 | -4.0711197 | 1.1094752 | 0.7056338 | 0.0987973 |
| ALB | -2.2266075 | -0.1452531 | -0.0732705 | -16.119826 | 0.8766075 | -2.3125831 | 0.7056338 | 0.6061826 | 0.1649558 |
| A.G | -0.2056375 | -0.0472918 | -0.0248308 | -9.825744 | 0.0663326 | 0.2007528 | 0.0987973 | 0.1649558 | 0.0825138 |
corrplot::corrplot(par_cor_1e$estimate, type = "lower", method = "color",
addgrid.col = "gray50", addCoef.col = "black",
title = "Partial correlation between all numeric variables one-vs-one.")
Figura S1.1. Mapa de calor con las correlaciones parciales entre todas las variables por pares. Las correlaciones parciales se han obtenido con la función ‘pcor()’ del paquete ‘ppcor’ y el gráfico se ha hecho con la función ‘corrplot’ del paquete ‘corrplot’.
pca <- princomp(x = data_num, cor = T)
fviz_pca_var(pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, axes = c(1,2), title="PC1 - PC2", labelsize = 4) +
fviz_pca_var(pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, axes = c(1,3), title="PC1 - PC3", labelsize = 4) +
fviz_pca_var(pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, axes = c(2,3), title="PC2 - PC3", labelsize = 4)
Figura S1.2.
library(scatterplot3d)
scatterplot3d(pc$scores[,1], pc$scores[,2], pc$scores[,3],
xlab="PC1", ylab="PC2", zlab="PC3",
color = ifelse(data_filt$Selector == "1", "green", "yellow"),
pch = ifelse(data_filt$Gender == "Male", 1, 15),
main = "3D PCA with the first 3 components colored by selector and gender.")
Figura S1.3. PCA en 3D con las primeras tres componentes principales, con los individuos coloreados por la variable selector y con distintas formas para la variable Gender.
data_filt %>%
dplyr::mutate(Selector = Selector %>% dplyr::recode("1" = "Patients",
"2" = "Controls")) %>%
melt() %>%
ggplot(aes(x=variable, y=value, fill=variable)) + geom_boxplot(outlier.colour = "Gray50") +
coord_cartesian(ylim = c(0,150)) +
facet_wrap(~Selector) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 0.95),
legend.position = "none") +
ggtitle("Variables of our study separated in patients and control subjects.") +
ylab("") + xlab("")
Figura S2.1. Gráficos de cajas de las variables numéricas presentes en los datos del ejercicio 1, separados en dos paneles según su condición (‘Selector’) de pacienteso controles. El eje Y ha sido limitado para una mejor visualización de las cajas a cambio de no visualizar algunos outliers de algunas variables.
Everitt and Hothorn (2011). An introduction to applied multivariate analysis with R. Springer.
C.M. Cuadras (2018). Nuevos métodos de análisis multivariante. CMC Editors.
Elsieo Martínez H. Tratamiento matricial de los datos multivariantes. Disponible online.
Kassambara (2017). Practical guide to principal component analysis in R. STHDA. Disponible online.
Pryiank Goyal (2020). PCA. Rpubs. Disponible online.